Monte Carlo simulations of two-dimensional charged bosons 
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Quantum Monte Carlo methods are used to calculate various ground state properties of charged 
bosons in two dimensions, throughout the whole density range where the fluid phase is stable. 
Wigner crystallization is predicted at r s ~ 60. Results for the ground state energy and the momen- 
tum distribution are summarized in analytic interpolation formulas embodying known asymptotic 
behaviors. Near freezing, the condensate fraction is less than 1%. The static structure factor S(k) 
and susceptibility x(^) ar e obtained from the density-density correlation function in imaginary time, 
F(k, t). An estimate of the energy of elementary excitations, given in terms of an upper bound 
involving S(k) and is compared with the result obtained via analytic continuation from F(k, r). 

PACS numbers: 02.70.Ss, 05.30.Jp 



I. INTRODUCTION 

The two-dimensional fluid of point-like spinless bosons 
interacting with a 1/r potential has drawn attention in 
the literature** as a model in quantum statistical mechan- 
ics which parallels the physically more relevant fluid of 
electrons. At zero temperature, the model is specified by 
the coupling parameter r s = 1/ 0ma b , where n is the 
density and the Bohr radius. For small r s the sys- 
tem is a weakly coupled fluid, well described by the Ran- 
dom Phase Approximation^ whereas it becomes strongly 
correlated and eventually undergoes Wigner crystalliza- 
tion upon increasing r s . Several results for the ground 
state energy, static structure, screening properties and 
elementary excitations have been reported using the Cor- 
related Basis Function theory^ 1 * 3 " various implementation 
of the Singwi-Tosi-Land-Sjolander (STLS) formalism, 4,5 
and the Overhauser modelA The momentum distribu- 
tion has been calculated for low r s in the Bogolubov 
approximation^ A comparison between the STLS results 
for the 1/r potential and the ln(r) potential has been re- 
ported by Moudgil et aim- 

Although the charged boson model may find applica- 
tions to superconductors, either as a system of bound 
electron pairs** or in terms of an effective action with 
Fermionic degrees of freedom integrated out^* no di- 
rect realization of the system is experimentally avail- 
able. Therefore numerical results provided by quantum 
Monte Carlo (QMC) simulations constitute the only reli- 
able benchmark for analytic approaches. Extensive sim- 
ulation results are available for 3D charged bosons*)***** 
and for the 2D system with the ln(r) interaction! 1 ^ 4 *^ 

In this work we present QMC results for several ground 
state properties of the 2D fluid of charged bosons with the 
1/r potential. We use two different algorithms, namely 
diffusion Monte Carlo (DMC)j»£ which is more efficient 
in the calculation of mixed averages, and reptation quan- 
tum Monte Carlo (RQMC)f»£ which gives easier access to 
correlations in imaginary time. The exact ground state 
energy and the mixed estimate**^ of the one-body den- 



sity matrix are calculated with the former. Unbiased 
estimates of the static structure factor and the suscepti- 
bility are instead obtained, using RQMC, from the auto- 
correlation in imaginary time of the density fluctuation 
operator. The inverse Laplace transform of the same 
auto-correlation function yields valuable information on 
the spectrum of elementary excitations. 

II. METHOD 

Quantum Monte Carlo is the method of choice for 
strongly interacting bosonic systems in their ground 
state, because it yields exact numerical results for a num- 
ber of quantities, subject only to known statistical errors. 

The DMC method 1 - samples a probability distribu- 
tion proportional to the "mixed distribution" f(R) = 
$(R)$>(R), where R = {ri, • • • ,r N } is a point in the 2N- 
dimensional configuration space of the system, ^f(R) is a 
trial wave-function, and is the ground-state wave- 

function. The exact ground state energy is obtained as 
the average over the mixed distribution of the local en- 
ergy, El(R) — 1 &(R)~ 1 H'$>(R). For a general operator 
not commuting with the Hamiltonian, ground-state aver- 
ages can be approximated by the extrapolated estimate 
(twice the average over the mixed distribution minus the 
variational estimate), 16 which leads to an error quadratic 
in the difference ($ — ^>). Our results for the one-body 
density matrix are given in terms of this extrapolated 
estimate, as in Ref. Il2l 

For operators diagonal in R we avoid mixed esti- 
mates resorting to the RQMC method 17 (one could al- 
ternatively use the forward walking technique 18 within 
the DMC method). In RQMC, the evolution in imag- 
inary time of the system is represented by a time- 
discrctizcd path X ~{Rq, • • • , Rm}- The algorithm sam- 
ples the distribution P(X) = H{R ) 2 Il^i 1 G{R i - 1 
Rf,e), where G{R — > i?';e) is a short-time approx- 
imation to the importance-sampled Green's function 
*(i?')(-R'|exp(-eff)|i?)4'(i?)- 1 . Assuming M is large 
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enough, the inner time slices of the path are individ- 
ually sampled from the distribution <£>(i?) 2 , and se- 
quentially sampled according to the quantum dynami- 
cal fluctuations in the ground state. Pure estimators, 
($|0|$) = ((O(Ri))), and imaginary-time correlation 
functions, c(r) = ((0(Ri)0(Ri+ n ))), are thus readily ac- 
cessible (here ((•)) means average over the random walk 
in the space of quantum paths X, and r = ne)H 

In all simulations we consider a system of N parti- 
cles in a square cell with periodic boundary conditions. 
The trial function is chosen of the pair product form, 
= exp(-^2ij(u\ri - r,!)), where u(r) is the RPA 
pseudopotential following Ref. [PJ Both the pseudopo- 
tential and the Coulomb interaction are evaluated us- 
ing generalized Ewald sums^ As usual jiiiiS we estimate 
the finite-size effect on the ground-state energy from 
variational Monte Carlo simulations. Variational ener- 
gies En, calculated with N in the range 25-200, are 
used to determine the best-fit parameter in the form 
Eoo = EN+a(r s )/N+b(r s )/N 2 . Assuming that the same 
size dependence holds for the exact DMC energies, the 
optimal parameters a(r s ) and b(r s ) are then used to ex- 
trapolate to the thermodynamic limit the result of a sin- 
gle DMC simulation with N = 52. Other quantities have 
comparatively smaller finite-size errors, typically below 
the statistical accuracy of the present simulations. 



III. RESULTS 

A. Ground-state energy 

The DMC ground state energies of the 2D bosonic 
fluid in the thermodynamic limit are compared in Table 
[I] with the results obtained with the Singwi-Tosi-Land- 
Sjolander (STLS) method by Gold^ with a parametrized 
wave function approach by Sim, Tao and WuA and within 
the Hypernetted Chain Approximation (HNC) by Apaja 
et al.. 1 While all computations agree qualitatively, we 
note that the agreement between HCN and the exact 
DMC results is particularly good. Our DMC results can 
be accurately reproduced by the parametrized function: 



E g {r s ) = ~[a r b s ° + + a 2 r b s 2 + 



a3r s 
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(1) 



where ciq and bo are fixed by the small r s behavior- 
(E(r s — > 0) ~ — 1.29355/r 2 ^ 3 ), b\ is fixed requiring a 

constant sub-leading term for r s — > 0, bi and 63 by re- 

1 3/2 

quiring leading terms in r~ and r s for r s — > 00. The 
final values of the parameters are c = 7/40, do = 0.2297, 
01 = 0.161, a 2 = 0.0594, a 3 = 0.01017, b = 80/21, 6 a = 
94/21, 6 2 = 73/14 and b 3 = 40/7. The reduced x 2 for the 
fit with 4 parameters and 7 data points is 1.5 at r s = 1. 
The above interpolation formula allows to obtain, by 
means of the virial theorem, the unbiased estimator of the 
average kinetic energy (ke) = —d(r s E g )/dr s as well as of 

the inverse compressibility 1/pKr = 
both reported in Table |U 
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FIG. 1: Ground-state energy for 2D triangular Wigner crys- 
tal (WC), bosons (B), unpolarized (UP) and polarized (P) 
fermions as a function of r s . Wigner crystal and fermion 
data are from Ref. I2H . On purpose of clarity we plotted 
rV 2 {E(r s ) — ci/r s )), with c\ = —2.2122, while the inset shows 
the corresponding E(r s ) curves. Points with error bars are 
size-extrapolated DMC results, continuous curves are analyt- 
ical fits. 



In Fig. Q our results are compared with the previ- 
ous DMC results by Rapisarda and SenatoreSS for 2D 
fermions and for the 2D Wigner crystal. In two di- 
mensions bosons crystallize at r s ~ 60 and fermions at 
r s ~ 34. The difference in critical density is analogous 
to the difference obtained in the 3D case, where bosons 
crystallize at r s = 160 and fermions at r s — 100>ii 



B. Momentum distribution 

The one-body density matrix n(r) and its Fourier 
transform, the momentum distribution n(k), have been 
computed performing random displacements of particles 
on the sampled configurations as explained in Ref. Il3l 

At variance with the 3D case^ the standard procedure 
leads to strong size effects due to the slow convergence 
of n(r) to its asymptotic limit no = lim r ^ oc n(r). We 
removed the size-effect adopting the correction proposed 
by Magro and CeperleyA^ for 2D bosons with In r inter- 
actions. Our results for the one-body density matrix are 
shown in Fig. (J2J. 

Extending to the 2D case the discussion presented for 
3D charged bosons in Ref. 0, we fix the divergence of 
the momentum distribution at small k 



n{k -> 0) 



n no^rs/2 



4S(k) (fcr ) 3 / 2 



(2) 



where n is condensate fraction, and r — r s a,B. The 
cusp condition 22 instead gives information on the short- 
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TABLE I: Ground state energy for bosons from VMC and DMC, extrapolated to the bulk limit and compared with estimates 
from approximate theories. We also give the average kinetic energy and inverse compressibility obtained from Eq. All 
values are in Rydberg per particle, the digits in parenthesis represent the error bar in the last digit. 



fa 


E (DMC) 


E (VMC) 


HNC 1 


STW 3 


STLS 5 




1/uKt 


1 


-1.1448(5) 


-1.14269(7) 


-1.1458 


-1.1062 




0.2903 


-0.531 


2 


-0.6740(2) 


-0.67192(6) 


-0.6740 


-0.6631 


-0.6484 


0.1442 


-0.3582 


5 


-0.31903(5) 


-0.317456(6) 


-0.3185 


-0.3133 


-0.3078 


0.04896 


-0.187 


10 


-0.17480(5) 


-0.17385(3) 


-0.1741 


-0.16685 


-0.1724 


0.01961 


-0.1097 


20 


-0.093387(8) 


-0.092903(3) 


-0.0928 


-0.086024 


-0.0959 


0.007533 


-0.06177 


40 


-0.048986(8) 


-0.048737(2) 








0.00286 


-0.03359 


75 


-0.026965(6) 


-0.0268246(8) 








0.001189 


-0.01892 
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FIG. 2: One-body density matrix n(r) at r s =1, 2, 5, 10, 20, 
40 and 75 



TABLE II: Best fit parameters for equation (QJ. The last line 
reports the value of g(0) from Fig. Q as used in the fit of 
n(k). 





1 


2 


5 


10 


20 


40 


75 


no 


0.531 


0.38 


0.176 


0.0677 


0.018 


0.001 


0.0007 


do 


0.839 


0.853 


0.475 


0.977 


0.861 


1.21 


2.59 


ai 


44 


3.5 


5.46 










a z 


-0.086 


0.492 


2.17 


1.96 


1.05 


0.946 


0.098 


03 


0.696 


0.56 


-2.1 


-1.13 


-0.08 


-0.74 


0.627 


tin 


1.13 


0.226 


0.23 


-0.01 


-0.163 


0.184 


-0.103 


a 5 


0.135 


0.192 


0.28 


0.7 


-0.006 


-0.014 


-0.024 


(Mi 


-111 


-6.07 


1.12 


-1.52 


0.849 


3.44 


0.576 


a 7 


6.98 


2.29 


1.45 


1.86 


2.61 


1.99 


2.59 


9Q 


0.21 


0.078 


0.01 
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no— 


0.537 


0.398 


0.230 











range behavior of the momentum distribution: 



n(k — > 00) 



4^ 2 .g(0) 
(fcr ) 6 



(3) 



where g(0) is the pair correlation function at r = 0. 
Moreover, at small density, we expect the momentum 
distribution to be approximately Gaussian, in agreement 
with harmonic theory for the crystalline phase. 

We have collected all this information in a fitting for- 
mula to interpolate the DMC data for the momentum 
distribution n(k): 



n(k) 



(2TT) 2 pn S 2 (k) 

0,2 



n ov /r s /2 



,3/2 



2 / 2 

-K /a Q 



+ a 3 + a 4 v k + as/* e 



where K — kr^. Given the known values of the density 
and of g(0) (see next section), we determined the remain- 
ing parameters by a least-squares fit to the DMC data 
on n(k), n(r) and on the average kinetic energy. 



Table HH contains the best-fit parameters and the re- 
sulting value of the condensate fraction rig . The conden- 
sate fraction decreases very rapidly with increasing r s , 
the depletion being already 50% at r s — 1, in agreement 
with the result of the Bogolubov theory^ (in 3Di£ a sim- 
ilar depletion occurs at r s = 5). For large couplings, the 
Bogolubov theory overestimates the condensate fraction. 
In a wide density range in the liquid phase, say r s > 20, 
no is of the order of 1% or less. Such small values, ob- 
tained by fitting Eq.|3]to the extrapolated estimates from 
the simulation, are presumably meaningful only as an in- 
dication of the order of magnitude. 



C. Imaginary-time correlation functions: static 
response function and static structure factor 

Information on charge response properties of the sys- 
tem like screening, plasma oscillations or polarization are 
contained in the imaginary time density-density correla- 
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FIG. 3: Static structure factor S(k) as a function of kro for 
r s — 1, 2, 5, 10, 20, 40, 60. Lines are only guide to the eyes. 



FIG. 4: The pair-distribution function for r s = 
1,2,5,10,20,40 and 60 (cubic spline interpolation of Monte 
Carlo data). Higher peaks correspond to higher values of r 3 . 



tion function: 



1 

2> 



e- TU S(k,u), (4) 



where Pk(r) = e'' k ' r ' and S(k,u)) is the dynamical 
response function. The correlation functions F(k, r) have 
been computed with RQMC for systems of 56 particles. 

The static structure factor S(k) is readily obtained 
from the imaginary-time density-density correlation 
function as: 



S(k) 



duj S(k,w) = F(k,0). 



(5) 



In Fig. J3J) we report the behavior of S(k) for various den- 
sities. As r s increases and approaches the crystallization 
density a sharp peak develops in correspondence with 
the first lattice wave-vector of the 2D Wigner crystal, 
kr = (2tt % /3) 1 / 2 ~ 3.3. 

In Fig. we report the pair distribution function: 



r)>. 



(6) 



At low density g[r) develops a high peak and long-range 
oscillations typical of a system approaching localization. 
As the density increases the effective repulsion between 
particles decreases and overlapping between charges be- 
comes possible. The behavior of S(k) and g(r) is quali- 
tatively in agreement with the findings of Apaja et a/., 1 
but for both functions the Monte Carlo results show more 
pronounced effects of correlations at low densities. 

The static response function x(k) can be evaluated 
from the relation 



X(k) = -2 



S{k,u>) 



doj = -2 



/ F{k,T)dr. 
Jo 



(7) 



In Fig. JSJ) we report the static effective interac- 
tion Vk/e(k,0) where Vk is the Coulomb interaction and 
e(k, 0) = 1/[1 + Vkx(k, 0)] is the static dielectric function. 
At low k the effective interaction is given the compress- 
ibility sum rule, 



1 



fc™e(fc,0) pK T ' 



(8) 



while in the short-wavelength limit it behaves like the 
Coulomb interaction. The minimum of u^/e(fc,0) deep- 
ens and shifts to larger k upon increasing r s . We note 
that a negative dielectric function cannot be interpreted 
as a signal of instability of the bosonic fluid due to the 
presence of the rigid background. As in the case of the 
structural properties, in the large coupling regime the 
Monte Carlo data for the effective potential show more 
pronounced features than the results of Apaya et al.. 1 
This is shown, in terms of the static response function 
X (k), in Fig. ©. 



D. Excitation spectrum 

The elementary excitations spectrum of the density 
fluctuation is contained in the dynamic structure factor: 



S(k,oj) = ^|(n|p fc |0)| 2 ^- W „o). 



(9) 



We estimate the energy dispersion of the collective 
excitation by fitting the imaginary time dependence of 
F(k,r) with F(k,r) = A(k)e- Ul ^ T + Biky-" 2 ^ . 
This amounts to represent the dynamical structure fac- 
tor S(k,uj) as the sum of two delta functions. When a 
single mode has a dominating spectral weight, its disper- 
sion u>i(k), is reproduced reasonably wellfil regardless of 
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FIG. 5: Effective interaction for r s = 2, 5, 10, 20, 40 and 60 
(open symbols, Monte Carlo data; lines, cubic spline interpo- 
lations). Deeper minima correspond to lower densities. The 
solid dots at k = are the values of 1/pKr from Tab. Q 

0.5 i ■ , ■ , ■ 1 



with their corresponding upper-bounds, at different den- 
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FIG. 7: The excitation spectra for r s = 2, 5, 10, 20 (full cir- 
cles and open diamonds) are compared with their respective 
upper-bound u>™ tn (solid lines). Dashed curves corresponds 
to data from Ref. [l| for r s — 5 and r s — 20. Curves with 
deeper minimum corresponds to lower densities. 
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FIG. 6: The static response function x(k) at r s = 1 (solid 
dots) and r s = 10 (open diamonds). The solid lines are from 
Apaja et al.^ 



the representation chosen for the remaining part of the 
spectrum (a delta function at u>2(k) in this case). 

Moreover, combining our results for x(k) and S(k) we 
obtain, by means of a sum-rules approachjiSiSi a rigorous 
upper bound for the plasmon dispersion: 



< 



2 P S(k) 



(10) 



At low k a single mode exhausts the sum rule. In this 
case, the upper bound in Eq. IjlQI) becomes an equality 
and the strength of the excitation coincides with S(k). 

In Fig. [7| we show our results for the excitation ener- 
gies extracted directly from F(k,r) and compare them 




FIG. 8: Excitation spectrum near the rotonlike minimum for 
r s — 10, 20, 40, 60. Full circles and open diamonds, data from 
two- exponentials fit to F(k, r); solid lines, upper-bounds from 
Ea.tTUl 



sities. On increasing r s a roton-like mode, close to the 
first reciprocal lattice vector of the Wigner crystal, devel- 
ops and softens. The evolution of this minimum as the 
crystallization transition is approached is shown in more 
detail in Fig. JSJl. 

In conclusions, we have presented an extensive QMC 
study of ground-state properties of 2D charged bosons. 
The present results constitute a valuable benchmark for 
theoretical approaches, showing their range of validity. 
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